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1 . Definition  of  the  Problem 


Nowadays  there  *s  a big  and  still  increasing  demand  for  fast  and 
precise  techniques  of  aerial  photo  interpretation.  For  many  sciences, 
authorities,  and  individuals,  the  use  of  information  gained  by 
aerial  photo  interpretation  is  the  most  reliable  and  up  to  date  basis 
for  their  tasks.  The  extraction  of  object  data  (such  as  the  shape 
or  the  location  of  an  object)  from  aerial  photos  is  related  to  the 
interpretation  problem.  Considering  the  problem  of  map  production  a 
restricted  but  well  known  set  of  objects  has  to  be  recognized  in  and 
extracted  from  aeriaj.  photographs  for  compilation  into  a map.  The 
visual  recognition  and  manual  extraction  of  the  contours  of  the  major 
roads  from  one  aerial  photograph  with  a size  of  9'  by  9'  takes  approx 
2 to  4 hours. 

An  effort  has  been  made  to  automate  this  problem  as  far  as  possible 
in  order  to  speed  up  and  ease  the  extraction  0f  data  from  photographs 
The  problem  of  recognizing  and  classifying  is  left  to  the  human  while 
the  time  consuming  task  of  finding  the  exact  contour  lines  will  be 
accomplished  semiautomatically . The  research  task  to  be  described 
below  can  be  stated  as  follows:  Given  grey  level  aerial  photographs 
taken  from  high  altitudes  at  scales  from  1:20  000  to  1:  74  000,  de- 
velop automatic  methods  for  use  in  an  interactive  image  interpre- 
tation system,  which  extract  the  coordinates  of  roads  fast  and  pre- 
cisely under  interpreter  control. 
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2.  Automatic  versus  interactive  extraction  methods 

A review  of  several  well  known  algorithms  for  automatic  object 
extraction  and  object  recognition  (contour  filtering,  contour 
tracing,  object  filtering,  correlation  techniques,  multispectral 
object  recognition  etc.)  revealed  that  a fully  automatic  solu- 
tion of  the  problem  is  not  yet  possible  mainly  for  the  following 
reasons: 

- the  algorithms  are  not  flexible  and  adaptive  enough 

- the  results  are  not  acceptable  in  case  of  noise  and  dis- 
tortions 

- the  classification  power  is  not  satisfying  in  case  of  complex 
real  world  objects 

- most  algorithms  require  too  much  computing  time 


Interactively  controlled  semiautomatic  methods  seen  to  be  more 
promising.  Clearly,  the  amount  of  interactive  control  to  be  pro- 
vided by  the  operator  (e.g.  an  image  interpreter)  should  be  mini- 
mized. The  following  procedure  of  interactive  extraction  of  roads 
from  photographs  is  proposed: 

- digitize  the  photograph 

- display  the  digitized  photograph  on  an  image  display  with 
interactive  devices  (track  ball,  light  pen  etc.) 

- recognize  visually  a (part  of  a ) road 

- indicate  the  starting  point  of  the  road  on  the 
display  screen 

- call  the  algorithm  to  be  used  for  the  extraction  of  the 
road 

- set  parameters  (e.g.  thresholds)  for  the  algorithm 

- let  the  algorithm  work  and  display  the  results  immediately 

- decide  visually  whether  the  results  are  acceptable  or  not 

- interrupt  the  algorithm  if  neccessary  and  decide  on  new 
steps 

- assemble  the  results  and  perform  manual  or  automatic  post- 
processing. 
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The  integration  of  the  interpreter  with  his  outstanding  ability 
to  solve  the  classif icatory  part  of  the  problem  is  considered 
to  be  a big  advantage.  In  order  to  speed  up  the  computational 
part  of  the  problem,  the  automatic  methods  are  not  very  complex. 
Three  methods  have  been  developed  to  accommodate  the  differing 
appearrance  of  the  different  roads.  All  methods  will  be  des- 
cribed in  principle, while  the  third  one  will  be  explained 
completely  and  in  detail  (see  also  [1,  2,  3]  ) . 


3.  Three  approaches  to  the  extraction  of  roads  from  aerial 
photographs 

3.1  Incremental  method 


The  incremental  procedure  with  adaptive  parameters  is  well 
suited  for  clearly  visible  roads  with  little  distortions 
and  occlusions.  It  is  not  affected  by  local  changes  of  the 
grey  level  or  width  of  the  road,  nor  by  crossroads  or  curves, 
as  it  follows  the  course  of  the  road  step  by  step  checking 
the  grey  level  situation  at  each  step  and  adapting  its  para- 
meter values  continuously  to  the  values  actually  found. 

The  interpreter  has  to  mark  the  starting  point  for  the  method 
at  a clearly  visible  part  of  the  road.  He  may  then  set  initial 
values  for  the  parameter's  grey  level  tolerance,  width  of  the 
road,  length  of  a step  etc.  Standard  values  will  be  used  as 
defaults  for  these  parameters.  Now  the  algorithm  tries  to 
follow  the  road  from  the  start  automatically.  Suppose  that 
the  points  pg  and  p^  have  been  found  on  the  road  (see  fig.  la) . 
The  algoritnm  continues  as  follows: 

- connect  pg  and  p1  by  a straight  line 

- prepare  the  next  step  by  marking  a point  p2  on  this  line 
at  a distance  of  the  initial  length  of  a step  beyond  p1 
(see  fig.  la) 

- evaluate  the  coordinates  of  all  points  of  a section  of 
a circle  through  p2 
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compile  the  grey  levels  of  all  these  points  into  a grey 
level  diagram  (cross  section  of  the  grey  level  situation 
at  point  p2,  see  fig.  1b) 


a)  b) 

figure  1 


- analyse  the  grey  level  diagram  to  discover  features  of  the 
road:  dark  road  with  bright  suroundings  ("valley"  in  the 
diagram) , bright  road  with  dark  suroundings  ("mountain"  in 
diagram)  etc. 

- compare  the  candidate  features  of  the  road  with  the  actual 
width  tolerance  of  the  road.  In  case  of  compliance  take  the 
location  of  the  grey  level  maximum  (or  minimum  in  case  of 
dark  road)  as  the  next  point  p2  on  the  road 

- if  there  is  more  than  one  candidate  which  fulfills  the 
specifications  of  the  features  of  the  road  decide  for  one 
of  them.  Criteria  for  this  decision  are: 

i)  the  direction  towards  the  end  of  the  part  of  the  road, 
if  the  end  point  has  been  set 

ii)  the  persistance  to  go  on  straight 

iii)  always  turn  right,  but  store  the  coordinates  for  later 
processing  of  the  other  candidates 

- if  the  road  proves  to  be  straight,  increase  the  length  of 
the  steps  until  the  next  curve  will  be  reached  (see  fig.  2a 
and  b) . A curve  of  the  road  implies  candidate  features  of 
the  road  which  are  not  located  near  the  centre  of  the  grey 
level  diagram  (see  fig.  3a  and  b) . In  this  case  the  length 
of  the  steps  is  decreased  automatically. 

- if  no  candidate  features  of  the  road  are  found,  ask  for  help. 


The  points  of  the  road  which  have  been  accepted  by  the  algorithm 
are  displayed  immendiately  for  visual  control  by  the  interpreter 
Correcting  and  completing  of  the  extracted  roads  can  readily 
be  performed. 


3.2  Binarizing  method 


The  area  binarizing  and  interpolation  method  attempts  to  cope 
with  disturbed  and  partially  occluded  roads.  To  initiate 
the  procedure,  the  operator  has  to  set  the  following  para- 
meters (e.g.  with  a trackball) : 


■p— ....I..— - , ,u  
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- starting  point  of  the  section  of  the  road 

- end  point 

- width  of  the  regie  a to  be  processed 

These  parameters  define  a rectangle  which  should  enclose 
completely  the  selected  section  of  the  road  to  be  extracted. 
Only  pixels  of  this  rectangle  will  be  processed,  regardless 
of  the  size  of  the  image,  thereby  reducing  considerably  the 
amount  of  computing  time.  It  is  an  additional  advantage  that 
the  operator  can  avoid  ambiguities  by  choosing  an  adequate 
rectangle. 

In  the  first  step  the  grey  level  matrix  of  the  rectangle  is 
binarized.  A cross  section  through  the  starting  point  defines 

the  grey  level  G and  the  width  W of  the  road.  The  rectangle 

^ 1 + — 
is  binarized  by  using  two  thresholds  Gr  and  Gr,  where  G^  is 

the  arithmetic  average  of  both.  All  pixels,  the  grey  level 

of  which  lies  within  this  grey  level  slice,  are  assigned  t^e 

binary  value  1,  all  others  are  set  to  0.  Fig.  4 shows  an 

example  of  a binarized  section  of  an  image. 


After  binarization  the  picture  consists  of  several  areas  not 
connected  to  each  otlnr.  The  next  step  is  to  eliminate  those 
areas  definitely  not  belonging  to  the  road.  One  criterion, 
that  is  used  to  decide  this,  is  the  shape.  All  parts  of  the 
road  are  presumed  not  to  exceed  the  width  W plus  a certain 
tolerance  value.  A measure  for  the  width  of  an  area  is  given 
by  ROSENFELD'S  distance  transform  The  distance  transform 

assigns  a value  to  each  pixel  of  an  area  according  to  the 
distance  of  this  pixel  from  the  nearest  border  of  the  area. 

After  the  distance  transform  has  been  performed,  all  areas  whose 
width  exceeds  W are  deleted  (see  fig.  5).  Areas  that  are  thinner 
than  are  preserved  because  parts  of  the  road  might  be 
partially  hidden. 


Fig.  5:  After  deletion  of 
big  areas 
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After  deleting  all  inappropriately  shaped  areas  the  next  step 
will  decide  which  of  the  remaining  areas  will  be  termed  parts 
of  the  road,  and  will  interpolate  the  gaps  between  these  parts. 

This  can  be  done  beginning  at  the  start  and  the  end  point  and 
proceeding  towards  each  other.  If  a gap  is  reached,  all  re- 
maining parts  in  the  nearer  environment  are  considered  as  candida- 
tes for  the  most  probable  continuation  of  the  road.  The  probability 
of  a candidate  is  based  on  the  following  criteria: 

- the  elongation  of  the  candidate 

- the  length  of  the  gap  to  the  candidate 

- the  direction  of  the  gap  compared  to  the  global  direction 
of  the  road 

If  the  interpolation  fails  to  select  the  appropriate  candidate, 
the  interpreter  can  correct  the  decision.  Thus,  with  inter- 
active aids,  the  binarizing  method  will  be  able  to  extract 
curved  and  distorted  roads,  as  long  as  the  grey  level  along 
the  roads  do  not  vary  more  than  a certain  amount. 


3 . 3 Random  method 

This  method  will  extract  the  coordinates  of  a section  of  a 
road  even  if  it  is  partly  invisible,  by  recognizing  only  a 
few  points  of  the  clearly  visible  parts  of  the  road. 

In  the  first  step,  the  operator  will  define  the  area  of  interest 
interactively  which  should  contain  completely  the  section  of 
the  road  to  be  extracted.  The  width  of  the  rectangular  area  of 
interest  may  be  small,  medium,  or  large,  depending  on  the 
curvature  of  the  road  (see  fig.  6a). 

This  rectangular  area  is  divided  into  n slices  of  equal  size. 
Points  of  the  road  will  only  be  searched  on  the  n+1  sample 
lines  intersecting  the  rectangle.  The  first  line  containing  P 


s 


p 


can  be  used  to  supply  Information  about  the?  actual  grey  level 
and  the  width  of  the  road  if  those  parameters  are  not  defined 
interactively.  These  two  parameters  arc  not  adapted  while  the 
algorithm  proceeds,  'ho  grey  level  diagrams  of  the  inter- 
secting lines  will  be  searched  for  parts  which  most  likely 
belong  to  the  road.  To  facilitate  this  search,  all  points 
of  the  sample  lines  are  transformed  into  a normalized  coordi- 
nate system  (see  fig.  6b) . This  is  done  for  two  reasons  mainly: 

- to  minimize  access  to  mass  storage,  as  the  complete  picture 
matrix  does  not  fit  into  the  working  storage  of  the  computer 

- the  search  procedure  for  the  grey  level  diagrams  of  the 
sample  lines  is  simplified,  if  applied  to  points  of  a norma- 
lized coordinate  system 


Fig.  6:  Rectangle  containing  the  road  (a); 
transformed  slices  of  the  rectangle  (b) 

The  search  for  parts  of  the  sample  lines  which  may  belong  to 
the  road  concentrates  on  the  detection  of  certain  features  of 


the  road: 

- the  contrast  between  the  road  and  the  surounding  area  at 
both  sides  of  the  road 

- the  constant  grey  level  of  the  pixels  on  the  cross  section 
of  the  road 

- the  correspondence  of  the  features  of  the  roads  at  adjecent 
sample  lines 


Those  parts  of  the  sample  lines  which  are  definitely  recognized 


as  points  of  the  road,  are  transformed  back  into  the  original 
coordinate  system  of  the  image  and  define  a polygon  to  approxi- 
mate the  course  of  the  road. 

Note  that,  due  to  the  interpolation,  it  is  not  essential  to 


evaluate  all  visible  parts  of  the  road.  Thus,  the  procedure  is 
able  to  extract  roads  which  may  be  disturbed  by  intersecting 
roads,  or  occluded  by  bridges  or  trees  (see  fig.  7).  The  number 
of  slices  within  the  rectangle  depends  on  the  length  of  the 
section  of  the  road  and  on  the  curvature  of  the  road  as  well  as 
on  the  accuracy  required  (see  section  4 for  a more  detailed 
description  of  this  method) . 


Fig.  7:  Slices  of  successful  recognition 
of  the  road  (see  arrows) 


4.  Detailed  description  of  random  method 


When  we  started  to  develop  a method  which  should  be  able  to  find 
roads  in  aerial  images  in  spite  of  the  road  being  partially  hidden 
or  invisible  we  guessed  that  we  could  base  on  two  features  of  a 
special  road  [5^  . 

a)  the  grey  values  G are  within  a reasonable  tolerance  T 

b)  the  width  W of  the  road  is  nearly  constant 


To  overcome  faint  or  invisible  parts  of  roads  and  to  reduce  the 
set  of  data  to  be  processed,  the  problem  was  modified  to  finding 
parts  of  roads  within  rectangular  areas  of  interest  (see  also  3.2) 
A proper  choice  of  these  rectangles  can  guarantee  small  width  to 
length  ratios  resulting  in  small  areas  with  only  a limited  number 
of  pixels  to  be  processed  (fig.  8) . 


A systematic  search  and  processing  of  data  within  the  rectangle 
could  now  be  started  to  find  all  points  on  the  road  (fig.  9a) . 
Yet  a good  part  of  the  result  is  redundant  as  the  curvature  of 
the  road  within  the  area  of  interest  is  small  (no  turns)  by 
definition.  In  other  w^rds:  A smaller  number  of  vertices  could 
as  well  approximate  the  segment  of  the  road  appropriately  (fig. 


For  this  purpose  only  a limited  number  of  samples  is  taken  from 
the  rectangle  not  considering  the  course  of  the  road  between  two 
adjacent  samples. 


Experience  has  shown  that  it  is  sufficient  to  consider  sample 
lines  instead  of  2D-sample  areas.  Thus,  the  two-dimensional 
problem  to  find  segments  of  roads  in  a rectangular  area  can  be 
reduced  to  finding  points  of  the  road  on  one-dimensional  sample 
lines  ( f ig . 9c) . 


a)  systematic  and 
complete  search 
for  points  on 
the  road 


b)  search  for 
points  using 
only  samples 


c)  same  as  b but 
samples  are  one- 
dimens  ional 


figure  9:  area  of  interest  containing  a road  segment 
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Advantages  and  limits  if  this  sampling  method  are  obvious.  Small 

gaps  accidentally  located  between  two  (sample)  lines  are 

not  even  registered.  Sampling  of  a road  crossing  may  possibly  fail 


to  yield  a point  on  the  road  (fig. 
results . 


E?"- 


0) ,but  does  not  yield  incorrect 


figure  10:  points  found  on  a partially  hidden  road  using  constant 
width  and  grey  value  criterion 


Road  segments  touching  areas  with  nearly  the  same  grey  value  cannot 
be  recognized.  Patches  with  size  and  grey  value  similar  to  those 
of  the  road  located  nearby  the  road  may  yield  incorrect  results. 

Considering  this,  a first  series  of  attempts  was  made.  We  started 
from  the  center  of  the  first  sampling  line  where  the  road  should 
pass  using  a cone  as  area  of  expectation  for  the  next  point  on 
the  road.  This  procedure  was  refined  by  adjusting  the  grey  value 
to  be  expected  to  the  grey  value  of  the  actually  found  point  on 
the  road.  Such  deviations  of  the  grey  level  along  the  road 
could  be  compensated  within  a limited  range.  If  no  segment  of  road 
was  found  the  range  of  expectation  was  searched  for  grey  level 
intervals  greater  than  W.  If  this  search  was  successful  the  road 
was  supposed  to  being  tangent  to  an  area  of  the  same  grey  value 
as  the  road.  This  method  had  3 serious  drawbacks  which  are  illustrated 
in  figure  11 
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a)  deadlock  situation  b)  interfering  c)  abrupt  changes  of  grey 

parallel  values  along  the  road 

contour 

figure  11 

A deadlock  situation  may  occur  in  road  following  if  a short 
appendix  is  more  distinct  than  the  road  itself  (a).  Close-by 
parallel  contours  may  cause  the  procedure  to  leave  the  wanted 
road  (b) . Abrupt  changes  of  grey  values  along  the  road  cannot 
be  processed  and  the  procedure  may  stop  (c) . 

Another  approach  was  made  which  proved  to  be  more  promising. 

It  is  described  in  detail.  Fundamental  differences  to  the 
method  described  above  are: 

- the  grey  value  at  a point  is  only  of  inferior  importance 

- the  grey  value  profile  in  a line  is  examined 

- possible  candidates  are  selected  (this  could  be  done  randomly) 
i.e.  it  does  not  matter  whether  line  number  1 is  processed 
before  line  number  1'  or  v.v.) 

- the  number  of  candidates  is  reduced  to  one  per  sample  line 
considering  connectivity  and  smoothness  between  neighbouring 
lines . 
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The  first  step  is  to  define  the  rectangle  and  the  number  of 
samples  to  be  taken.  For  this  purpose  the  operator  has  to 
identify  a point  on  the  road  which  serves  as  a starting  point. 


This  could  be  done  via  light  pen  or  trackball.  The  length  and 
direction  of  the  rectangle  is  introduced  to  the  system  by  means 
of  a second  point  which  may  not  necessarily  be  part  of  the  road 
The  width  of  the  rectangle,  the  number  of  samples,  and  the 
estimated  width  of  the  road  (number  of  pixels)  are  entered  via 
function  keyboard  or  as  numbers  from  the  teletype  (fig.  12): 


enter  name  of  image  file 

enter  image  size  (e.g.  1024x1024  pixels) 

enter  picture  scale  (f.e.  1:44000) 


provide  width  of  rectangle 

" number  of  sample  lines 

" width  of  road 

" 2 points  on  rectangle 


done  once  per  image 
changed  if  necessary 
once  per  rectangle 


Due  to  the  limited  amount  of  data  to  be  processed  it  will  generally 

be  possible  to  keep  all  sample  lines  of  a rectangle  simultaneously 

in  core.  This  allows  fast  and  repeated  processing  of  the  original 

grey  value  data.  For  this  purpose  a normalized  matrix  is  computed 

with  the  number  of  rows  1^  corresponding  to  the  number  of  lines 

and  the  number  of  columns  k.  corresponding  to  the  number  of  pixels 

P 

along  a line  (fig  13). 


figure  13:  normalized  matrix  M containing  the  sample  lines 

As  data  transfer  from  mass  storage  to  core  takes  the  bulk  of 
processing  time  care  has  to  be  taken  to  access  each  image  line 
only  once.  Intersection  of  an  image  line  i (fig.  13)  may  yield 
from  1 up  to  k^  intersecting  points  with  a sample  line  1 de- 
pending on  the  slope  a of  the  rectangle 

a < 45°  1 point  of  intersection  per  line 

45  < a < go0  k (k<k^)  points  of  intersection 

a = 90°  k ^ points  of  intersection 

Each  image  line  i may  additionally  intersect  with  more  than 
one  sample  line.  The  transformation  algorithm  has  to  evaluate 
all  these  points  of  intersection  with  line  i and  store  them  in 
the  normalized  matrix  M.  All  further  processing  is  based  on 
data  contained  in  M. 
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Evaluation  of  points  on  the  road  segment  is  done  mainly  in 
2 steps: 

a)  processing  single  lines;  each  point  on  a line  is 
assigned  a probability  ph  to 

belong  to  the  road  (horizontal  probability) 

b)  points  with  high  probability  ph  are  examined  whether 
acceptable  neighbouring  points  exist  in  other 
lines  (vertical  probability  pv) 

In  each  line  a single  point  with  maximum  probability  p 
= f (ph,  pv)  is  retained. 


Evaluation  of  horizontal  probability 


Each  point  of  a horizontal  line  is  compared  with  its  neigh- 
bours in  this  line.  Depending  on  the  local  profile  of  the 
grey  level  diagram  ph  is  evaluated. 


ph  = d2*  C 


(D 


d2  is  the  second  order  difference  of  the  actual  point  with  up 
to  3 neighbours  to  the  right  and  to  the  left  (fig.  14).  No 
more  neighbours  are  considered  if  a point  of  return  in  the  grey 
value  diagramm  is  reached  (first  order  difference  changes  sign) . 
Thus  the  contrast  of  each  point  with  its  immediate  neighbouring 
points  in  both  directions  is  considered  while  more  distant 
points  are  ignored. 
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f.e.  in  fig.  14: 

d2k  = (gk~gk-2)_ (gk+3~gk)  = 2gk" (gk-2+gk+3) 
if 

sign(gk-gk-1}  = sign(gk-rgk-2}  + sign(gk-2'gk-3) 
and 

Sign  (gk+1-gk)  = sign(gk+2-gk+1 ) = sign <gk+3-gk+2) 


The  absolute  value  of  d2  is  high  on  edges  or  on  (positive  or 
negative)  peaks  of  limited  width  as  shown  in  figure  15.  This 
reflects  best  situations  in  the  test  materia] , where  the  pro- 
file of  roads  is  in  most  cases  a peak  containing  only  a few 
pixels  (1"5  depending  on  picture  scale  and  appearance  of  the 
road)  . 


figure  15:  examples  for  maximum  absolute  second  difference  d2 
and  factor  C considerung  similarity  of  grey  values 
of  neighbouring  points 
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To  take  into  account  a local  constancy  of  grey  values  a factor  C 


is  introduced.  This  factor  may  be  used  to  increase  probability  ph 
if  a neighbouring  point  exists  either  at  k-1  or  k+1  with 
a grey  value  similar  to  pk  (fig.  15). 


The  factor  C is  used  with  an  upper  and  lower  bound  (fig.  17). 
The  upper  bound  is  chosen  as  c = 2 in  order  to  get  d2'*C,|k_k, 
d2"*C" I . . „ (fig.  16).  The  lower  bound  C=1  is  chosen  to 
retain  high  probabilities  for  peaks  (fig.  16). 


A ' ~ ^ C*f(lclll) 

• I 

1 1 $ * ? ; J lA^I-Ull 

figure  17:  factor  C as  function  of  first  order  difference 
between  neighbouring  points 
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As  will  be  described  below  it  will  be  possible  to  make 
assumptions  about  the  grey  value  gas  of  the  road  in  a certain 
line.  This  value  might  be  used  to  increase  or  decrease  proba- 
bility ph  by  comparing  the  actual  grey  value  of  a point  in  a 
line  gac  with  gas.  A factor  f may  be  computed  which  reflects 
deviation  of  gac  from  gas  (fig.  18).  This  factor  may  then  be 
used  to  alter  ph: 

ph  <£=  ph  • f (2) 


figure  18:  factor  f^  as  a function  of  grey  value  deviation 
from  a reference  grey  value  gas 


Points  with  the  highest  values  ph  are  considered  to  be  the  best 
candidates  to  belong  to  the  road.  Adjacent  high  valued  points 
are  merged  to  a single  candidate.  To  reduce  amount  of  data  to 
be  processed  in  the  next  step  only  the  best  five  candidates 
in  each  line  are  retained.  The  number  of  five  points  has  been 

found  to  be  reasonable  by  experience.  It  has  turned  out  that  if  the 
five  best  candidates  from  (f.e.)  100  points  in  a line  do  not 
involve  a point  on  the  road,  the  sixth  one  will  be  the  right 
one  in  only  about  3 to  5 percent  of  the  cases.  This  is  not 
worthwile  keeping  the  additional  sixth  candidate.  The  candidates 
are  stored  together  with  their  probabilites  ph  for  each  line, 
such  that  if  50  samples  have  been  taken  50x5=250  candidates  are 
retained  for  further  processing  (fig.  19). 
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candidate  with  index  k in  line  1 
has  probability  ph  to  belong  to 
the  road.  There  are  two  points 
in  1 with  higher  probabilities 


figure  19:  tabel  of  the  5 maximum  horizontal  probabilities 
per  line 


Figure  20  shows  an  example  for  the  five  candidates  per  line  in 
the  right  part.  In  the  left  part  these  candidates  are  overlayed 
to  the  normalized  matrix  M.  Brightness  of  the  candidates  is  pro- 
portional to  ph. 
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The  polygons  are  each  valued  according  to  whether 

a)  some  candidate  exists  in  a line  or  not 

b)  a triplet  of  candidates  including  the  one  in  line  1 
is  nearly  col  linear 

Mere  existence  of  neighbour inq  candidates  yields  a component 
pv  i for  the  vertical  probability  according  to 

4 

pvl  . I (5-n)Pl.n  ,3, 

n-  1 

p^_n  is  the  probability  of  the  vertex  of  the  polygon  in 
line  (1-n).  It  can  bo  seen  that  existence  of  candidates 
noar  line  1 (actual  candidate)  with  high  probabilities  P^_n 
yields  a good  value  for  pvl.  Example  in  figure  22  yields 

pvl'  - 4 • pA_j  + 2 * P1-3+P1-4 
pvl"  - 4 * P1_1+2  * P1_3+P,,1_4 

If  f.e.  candidate  with  index  k"  in  1-4  has  already  had  a 
"better"  neighbourhood  than  the  one  with  index  k'  tP" i_4>P' i_4 1 » 
pvl"  will  become  greater  than  pvl'.  In  this  way  vertical 
connections  between  candidates  are  proliferated  from  "top"  to 
"bottom"  resulting  in  increasing  probabilities  along  a chain 
of  points. 

To  take  into  account  smoothness  of  the  road  local  collinearity 
is  considered  by  testing  all  possible  triplets  of  candidates 
in  the  vertical  neighbourhood  which  include  the  actual  candi- 
date. This  results  in  a second  component  pv2  where 

n-1  4 

Pv2  " 1 , vm’  (Pa-m+Pa-n)/2 

m-1  n-2 


(4) 


vm  is  a measure  for  the  vicinity  of  a candidate  in  line  a-m 
to  the  straight  line  defined  by  candidates  in  lines  a and  a-n 
(fig.  23).  It  is  defined  as 


| 

n 


figure 

23:  horizontal 

deviation  d„  of  a candidate 
m 

v = I d -d  I 

for 

Id  1 < d 

m ' max  m 1 

1 m 1 max 

and 

vm  " » 

for 

Id  1 > d 

1 m1  - max 

(5) 

dm  is  the  deviation  of  the  candidate  from  the  straight  line  and 
d is  the  maximum  deviation  tolerated  (f.e.  cl  =3)  . Triplets 

Ilia  X IuaX 

with  deviations  greater  or  equal  dmax  will  yield  no  contributions 
to  pv2. 


dm  is  computed  as 


d 


m 


(ka-i.-"+(n-"0'kal/n-Vn> 


where  . m «•  1 , n-1  and  n ■ 2,4 


(6) 


The  measure  v^  is  multiplied  with  the  average  of  the  probabilities 
fo  the  two  points  whose  probabilities  we  do  already  know  (see 
e.g.  (4)). 
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The  probability  pv2  may  than  consist  of  up  to  6 terms  according 
to  5 lines  considered  (fig.  24) . 
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figure  24:  The  six  triplets  which  are  tested  for 
collinearity  in  5 lines 

As  stated  above  pvl  and  pv2  are  computed  for  each  combination 
of  points  within  the  cone  and  with  assumptions  as  made  above. 

The  vertical  probability  of  the  actual  point  is  then 

pv  = pvl  + pv2  (7) 

The  maximum  value  of  pv  is  finally  used  and  added  to  the  horizon- 
tal probability  ph  of  the  point  yielding 

p = ph  + pv  (8) 

This  is  done  in  forward  direction  starting  from  the  first  line 
until. the  last  line  is  reached.  A road  configuration  as  shown 
in  figure  25a  would  result  in  candidates  with  varying  probability 
p as  shown  in  figure  25b.  To  suppress  appendices  symmetrically, 
p has  to  be  computed  once  again  from  bottom  to  top  starting  from 
the  last  line  resulting  in  probabilities  shown  in  figure  25c. 





darker 
parts  of 
the  road 
correspond 
to  higher 
probabili- 
ties 


a)  candidates  with 
high  values  ph 


b)  high  values  p c)  same  as  b)  but  additio- 

after  top  to  bottom  nal  bottom  to  top  pro- 
evaluation of  pv  cessing 


figure  25:  example  for  forward  and  backward  evaluation  of  pv 


In  order  to  limit  the  growing  of  p along  a smooth  chain  of  points 
and  to  weight  influence  of  ph,  pvl  and  pv2  some  constants  have 
been  introduced,  pvl  may  contribute  10  portions  (eq.  (3))  while 
pv2  may  contribute  6 portions  only  (eq.  (4) , fig.  24) . On  the 
other  hand  collinearity  of  neighbouring  points  is  a better  hint  for 
road  segments  than  the  mere  existence  of  these  points.  It  has 
turned  out  to  be  useful  to  increase  influence  of  collinearity  by 
50  % . pv  and  ph  have  been  given  equal  influence.  Therefore, 

pvl  pvl/ 10 

pv2  +—  pv2/6  *1.5 
pv  (pvl +pv2) /2 

p (pv+ph)  /2 

results  in 

p=  ( (pv1/lO+pv2/6-1.  5)/2+ph)/2  (9) 

Equation  (9)  yields  probabilities  which  are  generally  high  or 
low  enough  to  decide  whether  a candidate  in  a line  is  the  very 
one  which  we  have  been  looking  for.  At  this  point  in  each  line 
the  candidate  with  the  highest  probability  is  selected.  Now  we 
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have  got  a polygon  with  one  vertex  in  each  line  approximating 
the  road.  Results  from  test  images  are  good  enough  not  to 
implement  sophisticated  methods  to  automatically  cornet  the 
course  of  the  road.  Instead, correction  is  accomplished  in  two 
simple  steps 

a)  single  points  within  a smooth  environment  are  corrected 
automatically 

b)  other  bad  results  have  to  be  corrected  interactively 
at  the  screen 


Correction  of  points 


5 neighbouring  candidates,  as  shown  in  figure  26,  are  considered.  As 
three  points  define  a second  order  curve, one  can  be  quite  sure  that 
a parabolic  segment  of  road  exists  if  even  four  neighbouring 
points  are  situated  on  a parabola.  A fifth  point  between  can 
then  be  adjusted  to  the  parabola  (figure  26) . In  all  examples 
shown  in  the  chapter  "results",  only  this  simple  automatic  method 
to  correct  points  has  been  used  if  not  stated  otherwise.  If 
there  are  results  which  cannot  be  tolerated  or  corrected  auto- 
matically the  user  in  our  implemented  system  has  the  possibility 
to  delete  or  update  parts  of  the  road  segment  explicitely.  For 
this  purpose  the  resulting  polygon  is  displayed  twice  in  full 
brightness:  against  a dark  background,  and  secondly 
overlayed  to  the  matrix  M as  was  shown  in  figure  20.  The  user  can 
decide  in  which  of  the  two  representations  he  will  delete  or 
update  points. 


figure  26:  correction  of  single  points 
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Assumptions  on  grey  values 

It  has  been  stated  above  that  knowledge  about  the  grey  value 
of  the  road  segment  could  be  used  to  increase  horizontal  pro- 
bability ph.  As  it  is  very  difficult  to  get  exact  grey  values, 
various  assumptions  on  them  are  made.  On  the  basis  of  these 
assumptions  different  results  (polygons  approximating  the  road 
segment)  can  be  evaluated  and  compared  to  select  the  best  of 
them.  This  proceeding  is  possible  and  reasonable  as  mere  com- 
puting time  to  get  an  approximating  polygon  in  a rectangle  of 
fifty  lines  with  100  points  each,  is  only  about  1 to  2 seconds. 
This  is  true  although  all  programs  are  written  in  FORTRAN  and 
are  not  optimized  to  a high  degree. 

In  the  programming  system  there  exist  two  modes  by  which  assump- 
tions about  the  grey  value  can  be  made. 

a)  an  assumption  is  made  by  the  operator  (interpreter) . 

The  corresponding  polygon  is  evaluated. 

b)  all  five  possible  assumptions  are  made  automatically 
and  successively  by  the  system  and  all  polygons  are 
evaluated  and  compared. 

Possible  assumptions  about  grey  values  are  as  follows: 

1.  The  grey  value  gas  (fig.  18)  of  the  road  is  assumed  to  be 
the  highest  one  existing  in  each  line  (road  appearing  bright 
compared  to  background) 

2.  gas  is  the  lowest  grey  value  in  each  line  (dark  road) 

3.  gas  Is  the  lowest  or  the  highest  grey  value,  that  means 
gac  is  compared  with  the  brightest  and  the  darkest  point 

in  the  line.  The  better  conformity  is  used  to  compute  f (2). 
This  assumption  is  good  if  grey  value  changes  abruptly  due 
to  different  surfaces  of  roads. 

4.  gas  is  the  grey  value  of  the  best  candidate  in  the  first  line. 
As  the  starting  point  is  defined  interactively  it  is  possible 
to  hit  a point  close  to  the  road  with  the  cursor.  Therefore 
only  a very  restricted  range  ground  the  cursor  position  is 
searched  for  candidates.  Only  the  best  one  (instead  of  5 
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in  the  other  lines)  is  taken  as  the  candidate  in  the  first 
line.  This  assumption  for  gas  is  good  for  homogeneous  grey 
values  on  the  surface  of  the  road. 


5.  No  assumption  for  gas  is  made.  Factor  f is  taken  as  f =1 

y y 


If  various  assumptions  have  been  made, the  best  result  has  to 
be  found.  This  can  easily  be  done  by  adding  all  probabilities 
p of  the  candidates.  The  maximum  sum  is  assumed  to  correspond 
to  the  best  result  as  it  reflects  best  appearance  of  the  pro- 
file of  the  road  and  coherence  between  candidates. 

The  overall  procedure  is  shown  in  figure  27. 


display  aerial  image 
definition  of  parameters 

I 

redefinition  of  parameters 
display  the  rectangle 

I 

loading  of  elements  of  Matrix  M 
into  core 


assumption  on  grey  value 

a)  a special  assumption 

b)  try  all  posibilities 


evaluation  of  probability  ph  for  all  points 
in  a line  retain  the  5 best  ones 

evaluation  of  probabilities  p=f(ph,pv)  for 
all  candidates  in  a line 

evaluation  of  probability  p=f(p,pv)  for  all 
candidates  in  a line  with  cone  rotated  1 80° 

t 1 

automatic  correction  of  single  points 

* 

repetition  if  more  assumptions  are  to  be  tested 

finding  of  the  best  result 

i 

display  the  result  for  interactive  correction 

repetition  with  new  assumption  if  result  is  not 

satisfying 

transformation  of  candidates'  coordinates  back  to 
image  coordinate  system 

* 

display  candidates  in  image 


figure  27s  overall  procedure  for  road  extraction 


for  all 
lines 


for  all  lines 
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for  all  lines 
from  bottom 
to  top 
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5 . Results 

In  this  section  results  are  discussed  which  were  obtained 
using  the  random  method  to  extract  segments  of  roads.  Sections 
of  photographs  with  a scale  of  1:74000  were  used  as  test 
material.  These  images  were  displayed  on  a Comtal  display 
with  a resolution  of  512  x 512  pixels  (9  bit  spatial  resolution). 
The  processing  was  done  using  a digitized  version  of  the 
same  image  with  a better  11  bit  resolution  (2048  x 2048  pixels). 
Roads  had  an  average  width  of  3 to  4 pixels.  Where  the  appearance 
of  the  road  was  good,  even  the  9 bit  resolution  was  good  enough 
to  provide  1 pixel  which  could  be  distinguished  from  the 
surrounding  of  the  road. 

Figure  28  shows  a section  displayed  on  the  screen.  A rectangle 
was  defined  in  the  upper  left  of  the  image  surrounding  a bright 
road  together  with  a darker  object,  possibly  a ditch.  Details 
can  be  recognized  better  looking  at  the  display  of  the  norma- 
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Fig.  30:  candidates  with  gas  Fig.  31:  resulting  approxi- 

being  bright  mating  polygon 


Fig.  32:  resulting  polygon 
overlayed  to  matrix  M 


In  figure  30  five  candidates  per  line  are  displayed.  The  brightness 
of  each  point  is  proportional  to  its  horizontal  probability  ph.  It 
can  be  seen  that  even  ph  for  itself  yields  good  results.  The  final 
result  is  shown  in  figure  31  and  32.  After  computing  the  probability 
p the  best  point  in  each  line  was  selected  and  automatic  correction 
of  single  points  done.  No  interactive  correction  by  the  operator 
was  neccessary. 


Fig.  33:  candidate  with  gas 
being  dark 


Fig.  34:  resulting  approxi- 
mating polygon 


: 


Fig.  35:  resulting  polygon  Fig.  36:  result  of  bright  road 

overlayed  to  matrix  M and  dark  ditch 


Figure  33  shows  the  equivalent  of  figure  30  but  a dark  appearance 
of  the  object  was  assumed.  The  ditch  near  the  road  was  extracted 
this  way  (fig.  34,  35)  and  transformed  back  and  displayed  together 
with  the  road  extracted  before  (fig.  36) . 
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The  next  sequence  demonstrates  extraction  of  a road  segment 
which  contains  3 crossroads  and  a darker  part  where  the  road 
is  no  more  visible  even  in  the  matrix  M which  was  taken  from 
the  2048  x 2048  image  (fig.  37,  38). 


Fig.  37:  definition  rectangle  Fig.  38:  normalized  matrix  M 

50  lines,  100  points  per  line 


Fig.  39:  resulting  polygon 


Fig.  40:  polygon  overlayed 
to  matrix  M 
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Figures  39  and  40  show  that  three  incorrect  points  remain 
in  the  dark  portion  of  the  image.  This  is  due  to  the  fact  that 
no  candidate  at  all  could  be  found  in  this  part.  Two  incorrect 
points  at  the  lower  end  of  the  polygon  are  caused  by  exceeding 
the  maximum  slope  of  the  polygon.  The  correct  points  were  dis- 
carded because  they  could  not  be  reached  by  the  cone  used  to 
compute  the  vertical  probability  pv.  Thus  the  interpreter  is 
forced  to  chose  a greater  number  of  sample  lines  or  to  shorten 
the  rectangle.  In  this  case  50  lines  had  been  taken  from  the 
rectangle  which  was  about  1200  pixels  long,  that  means  only 
every  24th  line  has  been  taken.  Nevertheless  it  was  easy  to 
correct  the  resulting  polygon  interactively  (fig.  41,  42,  43). 
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The  series  of  figure 
is  not  quite  as  good 
image  instead  of  the 


44  to  figure  r’0  shows  that  a result  which 
could  be  obtained  processing  the  512x512 
one  with  the  better  resolution. 


Fig.  44:  five  candidates  Fig.  45:  same  as  fig.  44  but 

per  line  overlayed  to  M;  50  lines, 

27  points  per  line 


Fig.  46:  resulting  polygon 
before  correction 


Fig.  47:  same  as  46  but  over- 
layed to  M 


If 
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after  interactive  correction  layed  to  M 


Fig.  50:  result  in  original  image 

Figure  50  shows  that  two  points  have  not  been  corrected  although 
they  are  located  beside  the  road.  Regarding  figures  47  and  49 
it  can  be  seen  that  decision  could  not  be  made  correctly  basing 
upon  the  display  of  M. 


Figures  51  to  54  demonstrate  application  of  the  procedure 
to  another  type  of  object,  a dark  river.  As  the  object  appears 
clearly  in  the  512x512  image  this  image  can  be  taken  as  the 
basis  for  processing. 


Fig.  51s  part  of  a river  to  be  Fig.  52:  resulting  approxi 
extracted  from  a 512x512  image  mating  polygon 


mrr>2 


Fig.  53:  resulting  polygon  over-  Fig.  54:  result  in  original 
layed  to  matrix  M image 
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The  final  example,  a road  crossing  the  river,  was  processed 
using  similar  parameters.  The  results  show  that  both  objects 
can  be  extracted  correctly  from  the  9 bit  resolution  image 
without  any  interactive  correction  as  at  least  one  contrasting 
point  per  line  exists  (fig.  55  to  58) . 


Fig.  55:  road  crossing 
the  river 


Fig.  56:  resulting  approxi 
mating  polygon 
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Fig.  57:  resulting  polygon 
overlayed  to  matrix  M 


Fig.  58:  result  in  original 
image 


6.  Conclusion  and  discussion  of  further  work 


Results  in  chapter  5 have  shown  that  the  procedure  des- 
cribed is  very  suitable  to  extract  roads  with  a "normal" 
appearance  in  the  image.  That  means 

- in  the  greater  part  of  the  sample  lines  at  least  one  pixel 
must  exist  with  reasonable  contrast  to  its  neighbourhood 

- extraction  of  closeby  parallel  highways  needs  modification 
of  the  procedure  as  only  a single  track  is  extracted  at 
the  present  state 

- varying  width  of  the  road  will  be  no  problem,  but  will 
possibly  yield  the  edge  of  the  road  rather  than  the  centre 
line 

- the  rectangle  and  the  number  of  sampling  lines  have  to  be 
chosen  such  that  the  slope  of  the  road  in  the  normalized 
matrix  M will  not  exceed  a limit  defined  by  the  aperture 
of  the  cone.  This  is  neccessary  as  a horizontal  probability 
can  only  reasonably  be  computed  if  intersection  of  the  road 
and  sample  line  is  unmistakable  and  candidates  in  adjacent 
lines  can  be  reached  by  the  cone  for  evaluation  of  vertical 
probability. 

Regarding  this,  the  procedure  can  as  well  be  applied  to  objects 
with  similar  appearance  as  roads  (rivers  f.e.).  Fig.  59,60,61,62 
show  how  ambiguities  near  the  object  may  cause  incorrect  results. 
This  may  happen  especially  if  the  defining  rectangle  is  very 
wide.  This  could  have  been  avoided  if  two  slender  rectangles 
had  been  chosen  instead  of  a single  wide  one.  Another  approach 
could  be  to  start  the  procedure  with  a slender  small  rectangle 
to  evaluate  a small  segment  of  the  road  (fig.  63)  and  start  the 
procedure  again  by  appending  automatically  a new  rectangle. 


Fig.  59:  defining  rectangle  Fig.  60:  resulting  approxi- 

for  a wide  area  mating  polygon 
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Fig.  61:  resulting  polygon 
overlayed  to  matrix  M 


Fig.  62:  display  in  original 
image  shows  3 incorrect  points 


Fig.  63:  defining  small  rectangle  Fig.  64:  first  rectangle 
for  incremental  operating  mode  automatically  defined 


Fig.  65:  further  step  in  Fig.  66:  final  step  in 

incremental  mode  incremental  mode 


The  new  rectangle  uses  the  endpoint  of  the  small  segment  already 
found  as  its  starting  point.  The  endpoint  is  evaluated  with  help 
of  the  last  direction  the  segment  had  taken  (fig.  64). 
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We  have  made  first  attempts  to  implement  such  an  incremental 
operating  mode  of  the  random  method  and  have  found  results 
promising  (fig.  63  to  66  and  67  to  69) . 
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In  both  examples  the  end  of  the  procedure  was  announced  by  the 
operator.  In  fig.  66  no  road  segment  could  have  been  found  if 
the  new  rectangle  is  defined  by  a simple  method  as  described 
above.  In  fig.  69  no  part  of  the  road  was  visible  beyond  the 
last  points  found. 


Another  attempt  which  could  be  made  to  automate  road  extraction 
to  a higher  degree  would  be  to  find  some  of  the  points  on  the 
roads  automatically.  These  points  could  then  be  used  as  starting 
points  instead  of  interactively  defined  points  at  the  present 
state. 
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